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ABSTRACT 

We study the orbital properties of stars in four (published) simulations of thick disks formed 
by: i) accretion from disrupted satellites, ii) heating of a pre-existing thin disk by a minor 
merger, in) radial migration and iv) gas rich mergers. We find that the distribution of orbital 
eccentricities are predicted to be different for each model: a prominent peak at low eccentricity 
is expected for the heating, migration and gas-rich merging scenarios, while the eccentricity 
distribution is broader and shifted towards higher values for the accretion model. These dif- 
ferences can be traced back to whether the bulk of the stars in each case is formed in-situ or is 
accreted, and are robust to the peculiarities of each model. A simple test based on the eccen- 
tricity distribution of nearby thick disk stars may thus help elucidate the dominant formation 
mechanism of the Galactic thick disk. 



1 INTRODUCTION 



Several mechanisms have been proposed to ex plain the formation 
of thick disks in galaxies (see Maie wskilll993h . However it is still 
unclear by which of these mechanisms thick disks preferentially 
form. This is despite the fact that more th an 25 years have passe d 
since it was first detected in the Milky Way dCilmore & Reidll983h . 
and that it has been established that this component appears to be 
ubiquitous in late-type systems (e.g. lYoachim & Dalcantonll20o3 . 
and references therein). 

Amongst the scenarios proposed to explain the formation 
of the thick disks are the direct accretion of stars from dis- 
rupted satellites (e.g. lAbadi et"^ i2003b). the thic kening of a 
pre-existing thin disk through a minor merger (e.g. | Quinn et al.l 



119931 : IVillalobos & Helmill2008l : iKazantzid is et SiV}2Qm. the scat- 
tering ormi^rationofst^^ Schonrich & Binnev I 
l2009al : iRoskar et alj l2008l : Ischonrich &"Binnevl |2009b), and in- 
situ trigered star formation during/aft er gas-rich mergers (e.g. 
iBrook et aI.ll2005l : lBournaud et al.ll2007|) . 

Even though studies of external galaxies have been fundamen- 
tal to establish the statistical properties of thick disks, it is likely 
that only for the Galactic thick disk we will be able to unravel 
its evolutionary path. For example, measurements of the phase- 
space coordinates for nearby thick disk stars allows reconstruction 
of their orbits, which contain imprints of the dynamical history, 
while their chemical abundances encode information about their 
sites of origin. Time is ripe to delve into more detailed predictions 
for the above-mentioned scenarios, because these have reached a 
level of maturity and detail that they warrant and permit a nearly 
direct comparison to observations. 



In this Letter we investigate how the orbits of thick disk stars 
can be used to distinguish between the various formation channels. 
In particular, we focus on the predicted eccentricity distributions. 
We expect our findings to be applie d soon to samples of nearby 
thick disk stars from SEGUE ( Y annv et al.l2009l : ISmith et alj200g) 
and RAVE (Steinmetz et al. 2006, Breddels et al. submitted), and 
in the long term, to the Gaia dataset which will provide much more 
accurate information for much larger samples of stars spanning a 
wide range of distances from the Sun. The orbital eccentricity-test 
should help to elucidate the dominant mechanism by which the 
Galactic thick disk formed. 



2 NUMERICAL EXPERIMENTS 

We have gathered four existing numerical simulations of late-type 
galaxies that, having all developed a thick disk component, clearly 
differ in the dominant formation mechanism. These are: 



(i) accretion and disruption of satell ites llAbadi et al.ll2003lj). 

(ii) disk heating by a minor merger ( Villalobos & Helmi 200^, 

(iii ) radial migration via resonant scattering dRoskar et al.l 
l2008l) , 

(iv) in-s itu formation during/after a gas-rich merger 
( iBrooketaU 2004 .2005) . 



2.1 Thick disk formation models 

Because each of the simulations mentioned above have been al- 
ready introduced in the literature, here we will only review their 
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Table 1. Relevant parameters for each of the simulations. (1) model, (2) Virial mass, (3) Bulge/Spheroid mass (4) Disk mass (thin+thick), (5) Thin disk 
scale-height (from a double power law fit to the vertical mass profile), (6) Thick disk scale-height, (7) Radial scale length of the thin disk, (8) softening lengths, 
(9) reference to the original articles. For comparison, estimated values for the Milky Way have also been included. For the cosmological simulations, the virial 
masses are defined as the mass enclosed within the radius where the local density falls below A = 100 times the critical density of the universe. A Hubble 
constant Ha = 70Mpc~^km/s is assumed when necessary. Disks scale-heights have been determined by fitting double-exponential laws to the vertical 
stellar density profile in the cyHndrical shell 2 < R/Rd < 3. 
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main relevant features, referring the reader to the original papers 
for further details. Table [T] summarizes their key parameters. 



2.1.1 Accretion scenario 



lAbadi et all ( l2003bl) showed that within the ACDM paradi sm, the 
accretion of stars from disrupting satellites in approximately co- 
planar orbits may give rise to an old thick disk component that 
comprises about one-third the mass of the much youn ger thin disk. 

In our sample we include the galaxy presented in lAbadi et al.l 

(l2003h) . which formed in a cosmological N-body/SPH simulation. 
This object, with a virial mass of the order of that of the Milky Way, 
was selected from a low-resolution simulation of a large volume of 
the Universe; and later re-simulated with much higher resolution. 
In this high resolution run, the mass per baryonic particle is ~ 3 x 
10^ Mq. The final mass for the thick disk in this galaxy (derived 
via a dynamical decomposition) is 1.1 x 10^^ Mq. 



2.1.2 Heating scenario 

In this model, a thick disk is formed by the dynamical heating that is 
induced by a massive satellite merging with a primordial, rotation- 
ally su pported thin disk. This scenario has been explored re cently 
by e.g. lVillalobos & Helmil l l2008h : lKazantzidis et alj \200A who 
have shown that 5:1 mergers and with a wide range of orbital in- 
clinations generate thick disks whose properties are in reasonable 
agreement with observations. In such a model, the bulk of stars that 
end up in the thick disk origin ate from the primordial di sk rather 
than from the accreted satellite l l Villalobos & Helmill2008l) . 

In our analysis we include one of the numerical experiments 
presented in Villalobos & Helmi (2008). In these simulations, the 
mass ratio between the satellite and the host is 0.2 and its initial or- 
bit is prograde and inclined by 30° with respect to the host disk. The 
mass per stellar particle in the simulation is rUp = 1.2 x IO^Mq, 
and the thick disk has a final mass of 1.2x IO^^Mq. It is important 
to clarify that only a small fraction of the thin disk component is 
present at the end of the simulation (~ 15 — 20% the mass of the 
original disk). This implies that, for this remnant to be the thick 
disk of a late-type galaxy, a new thin disk should form later from 
the cooling of fresh gas. This will lead to structural changes in the 
thicker component, which are not considered here. Nevertheless, if 
the growth of the new disk is adiabatic, then many characteristics, 
and in particul ar the eccentriciti es, are not expected to be dramati- 
cally different l lVillalobosl2009l) . 



2.1.3 Radial Migration scenario 

Stars in the thin disk may be trapped onto resonant corotation with 
spiral arms, and may migrate inwards and outwards along the spi- 
ral waves approximately conserving their angular momenta (and 
henc e eccentricity) and without leading to significant heating in the 
disk dSellwood & Binnevll2002h . However, since the vertical veloc- 
ity dispersion of stellar disks correlates with their surface bright- 
ness (Kregel et al. 2005J, the radial migration of stars from the in- 
ner regions (kinematically hotter) will result in the formation of a 
thicker disk component. 

Although this process has not been formally proposed as a 
thick disk formation scenario, numerical simulations suggest that a 
modest thick component may be buil t. Therefore, we incl ude in our 
sample the simulation presented in iRoskar et alj ( l2008h ran with 
the goal of characterizing the migration that takes place in galactic 
disks. The simulation starts with a dark matter halo of ~ 10^^ Mq 
where 10% of this mass is in the form of a hot halo gas compo- 
nent. This gas is allowed to cool and form stars self-consistently, 
mimicking the quiescent growth of disk galaxies, over a period of 
10 Gyr. The initial mass resolution is 10^ A/q for the baryons, and 
stellar particles have on average masses 3 x 



2.1.4 Gas-rich merger scenario 

The last scenario we explore consists in the formation of a thick 
rotating component during an active epoch of gas ri ch mergers in 
the past history of a galaxy jSrook et al.ll2004l . 20051) . This forma- 
tion channel differs fundamentally from the accretion model be- 
cause the bulk of thick disk stars are born in-situ rather than be- 
ing accreted from satellites. In this sense, this scenario might show 
certain similarities with the heating model. However, the latter re- 
quires the existence of a thin disk at early times, z ^ 1, in contrast 
to the merger scenario where the stars are already born in a hotter 
component. 

Here we an alyze the simulated galaxy introduced in 
lBrooketal.] ( l2004b . It formed in a semi-cosmological N-body/SPH 
simulation that includes heating/cooling of gas, star formation, 
feedback and chemical enrichment. Its dark halo has a quiescent 
merger history after 2; ~ 2, and a final baryonic content of the 
galaxy is ~ 5 x IO^^Mq (see Table [TJ. The mass per baryonic 
particle is ~ 2 x IO^Mq. The mass of the thick disk in this galaxy 
is ~ 2.2 x 10^ M0, identified as old stars (8.5 < age < 10.5) with 
relatively high rotation velocity (Kot > 50 km/s). 
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Figure 1. Vertical density profile of stars for each of the scenarios discussed 
in Section 12.11 The best fit mass- weighted double-exponential profiles ai'e 
shown with black solid lines, and the individual contributions of the "thin" 
and "thick" components are indicated by the red dotted curves. Notice that 
there is no significant thin disk in the heating scenario, thus only the ver- 
tical profile for the thick component is present. For the accretion scenaiio 
kinematical cuts have been applied in order to avoid contamination from the 
stellar halo (see text for details) 



The scenarios described in this Section are capable of pro- 
ducing a rotationally supported hot component whose properties 
resemble the 'thick disks' in galaxies. However, the relative pre- 
ponderance of such thick components does vary from galaxy to 
galaxy in our simulations. This is illustrated in Figure[T] where we 
show with solid dots the vertical mass profiles for each case in a 
cylindrical shell 2 < R/ Rd < 3; which minimizes the contri- 
bution from bars and bulges. Additionally, in the accretion model, 
particles identified as "spheroid" in Abadi et al. (2003b) have been 
removed (see Section 3). The error bars correspond to the rms ob- 
tained from one hundred bootstrap re-samples of the data, and they 
are generally smaller than the dot's sizes. The black solid lines 
show the best-fit double exponential profile found for each galaxy, 
together with its decomposition into the "thin" and the "thick" disk 
contributions, in red dotted lines. This decomposition is generally 
robust, although correlations exist between the relative density of 
each component and the scale-height of the thick disk. The scale- 
heights 20 obtained by minimizing are quoted in each panel, and 
the typical errors are of order ~ 5 - 10% for z^'"'" and ~ 20 - 30% 
for the thick component. For the migration scenario, the uncertain- 
ties are larger due to the stronger dominance of the thin disk, and 
are 15% and 60% for Zq'*'" and zo, respectively. Nevertheless, in 
all cases, the results presented below are robust to changes in the 
value of Zo within the uncertainties. 

Differences in the relevance of the thick component depend 
not only on the net efficiency of the respective formation process, 
but may also be influenced by the different initial conditions and 
simulation techniques (e.g. only the accretion and merger scenar- 
ios actually account for the full cosmological framework). Because 
of these basic differences between simulations, global properties 



such as the mass, rotation and size of each formed thick disk are 
expected to be diverse. Our focus, however, is on contrasting the 
specific dynamical properties of the stars in the thick component 
for each case, and in particular, their orbital eccentricities, which 
as we shall see below, are fundamentally related to the physical 
mechanism by which this component was built. 

In order to facilitate comparisons between the thick disks in 
our galaxies and in particular, also to that of the Milky Way, we 
re-scale the radial and vertical distances of the stellar particles in 
each galaxy by their corresponding thin disk scale-lengths and thick 
disk scale-heightsQ. In what follows, we will focus our attention 
on "solar neighbourhood volumes", equivalent to cylindrical shells 
between two and three scale-radii of the thin disk (2 < R/ Rd < 3). 
For comparison, RQ/Rd ~ 2.2 — 2.4, assuming a scale radius of 
7?d = 3.5 kpc for the Milky Way. 

2.2 Modelling of the orbits 

Kinematical surveys such as RAVE, SEGUE and ultimately Gaia 
provide phase-space coordinates of stars around the position of the 
Sun. This instantaneous information may be used to recover their 
plausible past orbits. This requires modelling the (unknown) Galac- 
tic potential, and possibly its evolution, which implies that the or- 
bital parameters derived for each star generally suffer from a certain 
degree of uncertainty, even if measurement errors are neglected. 

On the other hand, our numerical simulations allow us to track 
in time each particle, with full orbits that are known. Neverthe- 
less we prefer to mimic observations, and therefore we use the 
present-day position and velocity of each stellar particle as initial 
conditions for the integration of their orbits in the best-fit poten- 
tial of their host galaxy. We model each galax y as a four compo- 
nent system with an NEW i Navarro et al.lll997[) dark halo, a Hern- 
quist profile i Hemquis^ I990l) for the bulge and two Miyamoto- 
Nagai disks jMivamoto & Nagailll975h corresponding to the thin 
and thick disks contributions. The mass associated with each of 
these components is known for each simulation (see Table[Tll, and 
their various scale-lengths are chosen by requiring a good match 
to the circular velocity profile of the system up to a distance of 
20 kpc0. This is along the lines of previous work, where the circu- 
lar velocity of t he Milky Way is o ften used to constrain the model 



lar velocity oi t he Milky Way is o ite 
parameters (e.g. lHelmi et al.l2006f) . 
We define eccentricity as e 



center distance of the last orbit of each particle. With this defini- 
tion, for a circular velocity curve modeled with ~ 10% accuracy, 
the eccentricities obtained by numerical integration show a scatter 
±0.1 — 0.2 around their true value with no systematic trends (see 
Fig.|2j. Larger deviations are found as the eccentricity increases. 



3 RESULTS 

Baryons in galaxies are generally sorted in several components: 
a bulge, a disk (thin-l-thick) and a more extended and diffuse 

^ For the simulation by Villalobos & Helmi (2008) we assume the scale- 
length to be that of the Milky Way thin disk: Rd = 3.5 kpc, although our 
conclusions do not fundamentally depend on this choice. 
^ For the migration and merger scenarios only the total (thin + thick) mass 
of the disk is known. For these cases the relative mass ratio between the thin 
and thick components is also a free parameter 
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Figure 2. Eccentricity distribution of all stellar parti cles in a cylindrica l 
shell with 2 < R/Rd < 3 for the accretion scenario iAbadi et alj2003bl) . 
The various panels are for different heights above/below the plane (normal- 
ized to the thick disk scale-height, za). The thick solid black line shows the 
total distributions per z-bin, while blue-empty and red-shaded histograms 
distinguish between in-situ and accreted stars. The effect introduced by the 
numerical integration of the orbits is rather small: the e-distribution ob- 
tained from direct tracking of the particles orbits in the simulation is shown 
in dotted black line. The fractional contributions from the thin disk (ma- 
genta dot-dashed), thick disk (green long-dashed) and the spheroid (soHd 
black) to each eccentricity bin are shown below each histogram. The num- 
ber of particles A' included in each box is also quoted. 

spheroidal distribution, the stellar halo, that might extend well be- 
yond the luminous edge of the disk. Although each of these com- 
ponents has defining characteristics, some of its properties may 
change smoothly from component to component. For example the 
eccentricity of the orbits changes as we move away from the disk 
plane into the realms of the stellar halo. This can be seen in Fig- 
ure |2l where the eccentricity distribution of all stellar particles 
within a cylindrical radii 2 < R/Rd < 3 is plotted for different 
heights above/below the plane. Vertical distances are normalized 
to the scale-height of the thick disk, zq. The eccentricity distribu- 
tions for stars formed in-situ (defined as those bom within a dis- 
tance of 20 kpc from the main progenitor) and for those accreted 
are given by empty long-dashed blue and shaded red histograms re- 
spectively. For comparison, the eccentricity distribution measured 
directly from the simulation (i.e. by tracking individual particle or- 
bits) is given by the dotted histogram. This shows that no system- 
atic errors are introduced by the orbital integration in the model 
host potential. 

The lowermost bin \z/zq\ < 0.5 is largely dominated by thin 
disk stars with circular motions, as can be seen by the strong peak 
around e 0. 15 in the top left panel of Figure[2] These stellar parti- 
cles formed in-situ through local conversion of gas settled in a disk, 
into stars (blue empty histogram) within the main galaxy. As we 
move away from the plane the thick disk gains importance and the 
eccentricity distributions become flatter as the fraction of accreted 
stars increases. Further above and below the plane, the distribution 
is dominated by particles from the spheroidal component, with even 



higher characteristic eccentricities, i.e e ~ 0.7. These changes in 
the relative preponderance of the thin, thick and spheroidal com- 
ponents can be seen as their fractional contribution (magenta dot- 
dashed, green long-dashed and black solid, respectively) to each 
eccentricity on the bottom panel of all |z|-b ins. Here, we have u sed 
the dynamical decomposition performed in lAbadi et alj ( l2003bl) to 
assign stars to a given component. 

This galaxy, introduced in Abadi et al. (2003a,b), has a large 
stellar spheroid that contains more than ~70% of the total luminous 
mass. Figure |2] shows that its contribution dominates the high ec- 
centricity bins at all heights above/below the plane. To highlight the 
properties of the thick disk, and also to avoid confusion on the in- 
terpretation of our results and its comparison with other simulated 
galaxies (lacking such a prominent spheroidal component), we will 
exclude in what follows any stellar particle that have been assigned 
to the spheroid by the analysis performed in Abadi et al. (2003b). 

Figure[3]shows how the eccentricity distributions vary accord- 
ing to the different formation channels of the thick disk. Each panel 
corresponds to one particular model: accretion (top left), heat- 
ing (top right), migration (bottom left) and merger (bottom right). 
When relevant, the contributions from stars formed "in-situ' or 
''accreted" have been highlighted. In order to minimize the con- 
tribution from thin disk stars, we have focused on the vertical bin 
1 < \z/zo\ < 3. On the other hand, to avoid contamination from 
the spheroids in our simulations (this is unlikely to be important 
for the Galactic stellar halo because of its very low density), we 
only consider stars with rotational velocity > 50 km s^^. This 
corresponds to average azimuthal velocity at which there is a clear 
excess of stars with velocities larger than this threshold compared 
to the distribution at < —50 km s^^ in our simulations. Un- 
der the assumption of a non-rotating spheroidal component, this 
criterion will minimize the contribution of the stellar halo in our 
samples. Nonetheless, we have checked that our results are not 
strongly sensitive to this assumption. The v^f, > 50 km s^^ cut 
removes 33 per cent of the stars in Brook's model, but has a neg- 
ligible effect in the simulation by Villalobos & Helmi and Roskar 
et al. due to the suppressed contribution of satellite accretion. Re- 
call that for the Abadi's galaxy we have removed all the spheroid 
identified by the dynamical decomposition. Cuts on the cylindrical 
radii (2 < R/Rd < 3) also help to elude the contributions from 
central bars and bulges. 



Figure[3]shows that the eccentricity distributions of stellar par- 
ticles in these "solar neighbourhood" regions, and between one and 
three thick disk scale-heights above/below the plane, are different 
according to each model. For the accretion scenario, the distribu- 
tion is very broad, with a median eccentrici ty (e) ~ . 5 (in good 
agreement with the accreted component in iRead et alj|2008h . On 
the other hand, the heating of a pre-existing thin disk by a minor 
merger gives rise to a bimodal distribution. The dominant peak is 
at low eccentricity e ~ 0.2 — 0.3 and associated to the stars from 
the progenitor disk, while the second peak at e ~ 0.8, is brought 
by the disrupted satellite. Radial migration tends to preserve the 
initial (low) eccentricity distribution, with only one peak present at 
(e) ~ 0.2 and with a sharp cut-off at e ~ 0.6. Finally, in the merger 
scenario a prominent peak around e ~ 0.2 is visible which is as- 
sociated to stars formed in-situ during the epoch of active gas-rich 
mergers. But like for the accretion scenario, accreted stars from 
infalling satellites contribute to a high-eccentricity tail. 

Interestingly, Figure [3] shows, as expected, that stars formed 
in-situ have low eccentricity orbits regardless the mechanism (heat- 
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Figure 3. Comparison of the eccentricity distributions of each thick disk 
formation model for stars in the range 1-3 (thick-disk) scale-heights and 
cylindrical distance 2 < R/Rd < 3. The color and line coding are the 
same as introduced in Figure|2] 



ing, migration, merger) that places them at their current height 
above/below the plane. On the other hand, accreted stars from satel- 
lites always dominate the high eccentricity tails of the distributions. 
Changes in the relative proportion of in-situ versus accreted stars 
drive the differences seen in the histograms of each thick disk for- 
mation scenario. In other words, the analysis of the stellar eccen- 
tricities off the plane helps to unravel whether the bulk of stars was 
formed locally in the main progenitor (heating, migration, merger) 
or, on the contrary, was accreted from infalling satellites (accre- 
tion). 

Kolmogorov-Smirnov tests performed over randomly gener- 
ated subsamples of stars selected from each simulation show that 
^ 150 stars are enough to distinguish at 90% confidence level 
between all the scenarios investigated here. Moreover, in samples 
containing more than ~ 550 stars, the probability that the stars are 
drawn from the same distribution is found to be lower than 1%. 
Clearly a smaller number of stars (~ 200) is sufficient to distin- 
guish the accretion model from the rest. Note, however, that these 
estimates are only indicative since they are derived for particular 
realizations of a general class of models. 

It is important to recall that each of these simulations has pro- 
duced a galaxy with a different morphology. Furthermore the sim- 
ulations representative of the heating scenario as well as that of 
migration have, by construction, a suppressed contribution from 
accreted satellite galaxies. Nevertheless, the general behaviour of 
in-situ vs accreted populations are robust to the different simula- 
tions idiosyncrasies and can be traced back to the different phys- 
ical mechanisms related to where and how the stars were formed. 
Therefore, we do not expect the global properties of the eccentric- 
ity distribution (e.g. bimodality, high eccentricity tails associated 
to accreted populations) to differ significantly, but only in the de- 
tails, in other realizations of the same models. Although we have 
focused on "solar-neighbourhood" regions (2 < R/Rd < 3), our 
conclusions do not depend fundamentally on this choice. 



4 CONCLUSIONS 

In this Letter we have analyzed four numerical simulations of 
galaxies hosting a thick disk component of fundamentally differ- 
ent origin: (i) accretion of satellites, (ii) heating of a pre-existing 
disk by a 5:1 mass-ratio merger, (Hi) radial migration by resonant 
scattering and (iv) gas-rich mergers at high-redshift. 

We have compared the eccentricity distributions predicted by 
these different models for stellar particles in the "solar neighbour- 
hood", i.e. located in a cylindrical shell of radius 2 < R/Rd < 3 
and with heights 1 < \z/zo\ < 3. Thick disk stars formed in-situ 
have low orbital eccentricities e ~ 0.2 — 0.3, independently of 
the mechanism that brought them high above/below the plane: gas- 
rich mergers, heating or migration. On the other hand, and again 
regardless of the particular model, accreted stars always dominate 
the high-eccentricity tail of the distributions. Therefore, the char- 
acterization of the eccentricity distribution of the thick disk can be 
used to establish if this component was formed by the accretion of 
satellite s (Abadi et al. 2003) or alternatively locally w i thin the main 
galaxy (Br ook et al.ll2005Mvillalobos & HelmilllOOSi : iRoskar et al.l 
I2OO8) . However, given the various limitations of our set of sim- 
ulations, we cannot claim that it will be possible to make an un- 
equivocal classification among models (i) — (iv) based only on 
eccentricity. 

Nevertheless, the differences between the orbital eccentricities 
of in-situ and accreted stellar particles are encouraging in view of 
the various kinematic surveys mapping our Galaxy today and in the 
near-future. We believe that with a reasonable guess of the Milky 
Way potential, the analysis of the eccentricity distribution of thick 
disk stars at approximately 1-3 scale-heights should shed light on 
the formation path of the Galactic thick disk. 
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